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INTRODUCTION 

Chief,  Tedhsleal  XnfonaatlonOlvlsion 

The  Moving  finite  element  (MFE)  method  Is  a  new  approach  for  numerically 
solving  partial  differential  equation  (PDE)  systems;^'2*2)  It  Is  parti¬ 
cularly  well  suited  for  resolving  POE  solutions  which  may  contain  large,  mul¬ 
tiple  gradients  over  highly  disparate  scales  In  both  space  and  time.  These 
types  of  POE's  abound  In  such  basic  technical  disciplines  as  aerodynamics 
(with  emphasis  on  shear  layers,  shocks  and  their  possible  Interactions),  com¬ 
bustion,  plasma  physics.  Interface  phenomena,  continuum  mechanics,  and  other 
transport  processes. 


In  the  MFE  method,  grid  co-ordinates  themselves  are  dependent  variables 
which  are  calculated  continuously  at  each  time  step  In  order  to  minimize  PDE 
residuals.  This  feature  has  successfully  suppressed  numerical  dissipation  to 
very  low  levels  and  has  resolved  accurately  In  1-0  those  physical  transport 
processes  which  may  occur  over  extremely  small  scales  (e.g.  those  described 
by  the  Navler -Stokes  equations  for  viscous  compressible  fluids)  simultane¬ 
ously  with  other  physical  processes  which  may  span  macroscopic  scales.  The 
objectives  of  the  present  research  effort  are  to  extend  the  MFE  method  to  2-D 
and  to  Investigate  Its  basic  properties  and  needs  for  solving  Important  PDE's 
In  2-D.  The  work  tasks  which  are  contributing  to  the  achievement  of  these 
objectives  are  discussed  In  the  next  section  of  this  report. 


It  was  established  during  the  first  year  of  this  study  that  the  MFE 
method  does  extend  logically  and  practically  to  2-0.  An  Initial  experimental 
2-0  MFE  code  version  was  developed  In  that  early  work  for  the  purpose  of  con¬ 
ducting  continuing  theoretical  and  applied  mathematical  research  In  diverse 
scientific  contexts.  Also  In  the  first  year  of  study,  effective  MFE  node 
movement  properties  and  significant  node  savings  were  demonstrated  for  simple 
—  but  yet  significant  —  test  problems  in  2-0.  The  MFE  code  structure  was 
found  to  be  amenable  to  vectorlzatlon  and  to  use  on  envisioned  advanced  com¬ 
puters. 


Having  established  these  milestones  In  the  first  year,  work  during  this 
second  year  of  study  has  moved  to  a  new  level,  focussing  now  on  a  host  of 
those  issues  and  research  needs  which  are  essential  to  the  use  of  the  MFE 


method  In  truly  large-scale  2-D  computations.  As  In  the  past,  some  of  the 
emerging  MFE  results  In  2-D  may  be  unique  and/or  contrary  to  conventional 
wisdoms.  In  several  areas,  they  suggest  some  Important  lines  of  research 
which  should  perhaps  be  pursued  more  Intensively  In  the  future,  both  by  our¬ 
selves  and  by  others. 

This  report  Is  organized  as  follows:  the  Technical  Progress  section 
provides  a  substantive  statement  of  significant  accomplishments  and  progress 
toward  achieving  the  research  objectives.  The  Publications  section  provides 
a  emulative  chronological  list  of  written  publications  In  technical  journals 
and  Includes  those  In  press  and  manuscripts  of  articles  in  preparation.  The 
Personnel  section  lists  the  professional  personnel  associated  with  the  re¬ 
search  effort,  and  the  final  section  entitled  Interactions  lists  coupling 
activities. 

I.  TECHNICAL  PROGRESS 

Work  during  this  period  addressed,  among  other  tasks,  MFE  node  movement 
properties,  ODE  solvers  for  PDE  methods,  linear  systems  solvers  for  the  MFE 
method,  regularization  techniques  for  the  MFE  method,  boundary  conditions, 
alternative  co-ordinate  systems,  and  interface  tracking.  This  was  essentially 
a  building  year;  l.e.,  one  of  building,  testing,  and  Implementing  both  new 
and  existing  mathematical  methods  which  are  needed  In  the  Investigation  of 
the  MFE  method  In  2-D.  In  several  areas  above  our  MFE  work  has  barely  opened 
entire  new  lines  of  Inquiry,  and  results  are  either  preliminary  or  Incom¬ 
plete.  In  such  cases,  the  discussion  of  progress  Is  necessarily  general.  On 
the  other  hand  there  are  some  significant  Implications  which  are  already 
apparent  In  some  of  the  new  research  areas  which  have  been  broached  during 
this  period  for  which  the  potential  significance  Is  Indicated  more  specifi¬ 
cally  —  but  with  the  understanding  that  a  great  deal  of  additional  work  Is 
needed  In  order  to  further  substantiate  or  refute  the  early  Indications. 

1.  MFE  Node  Movement  Properties 

Burger's  test  example  cam  be  formulated  to  pose  computational  challenges 
of  varying  degrees  of  difficulty  to  PDE  methods  In  2-D.  This  can  be  done  by 
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selecting  Initial  conditions  which  give  rise  not  only  to  steep  shocks  with 
planar  profiles  (as  In  work  during  year  1  of  this  study)  but  also  to  highly 
skewed  wave  profiles  with  gradients  of  greater  of  lesser  magnitudes. 


The  PDE's  for  this  skewed  Burger's  model  problem  are  given  by: 

“t  *  -uux  -  vuy  ♦  v<uxx  ♦  “yy>  »> 

Vt  -  -UVX  -  VVy  +  V(VXX  +  Vyy)  (2) 

where  u  Is  the  x  component  of  velocity  and  v  Is  the  y  component,  and  v  Is  an 
effective  diffusion  coefficient.  Shocks  are  generated  with  gradient  magni¬ 
tudes  on  the  order  of  1/v.  Initial  conditions  which  produce  a  doubly  skewed 
wavefront  profile  are  shown  schematically  In  Figures  1  and  2.  (The  counter- 
posed  Initial  velocity  fields  are  designed  to  create  an  evolving  shock  pro¬ 
file  which  Is  skewed  In  both  the  x  and  y  components  of  velocity.)  Boundary 
conditions  are  given  by: 


u(0,y)  «  u(l,y)  -  0 
vx(0,  y)  *  vx(l,  y)  -  0. 
u(x,  1)  ■  0.2  sin  irx 
u(x,  0)  *  -.2  sin  irx 
v(x,  1)  ■  -1.  +  0.2  cos  irx 
v(x,  0)  ■  1.  +  0.2  cos  irx 


0.  <  y  <  1. 

0.  <  y  <  1. 

0.  <  x  <  1. 

0.  <  x  <  1. 

0.  <  x  <  1. 

0.  <  x  <  1  . 


The  MFE  nodes  are  fixed  by  zero  Dlrlchlet  conditions  along  the  top  and  bottom 
horizontal  edges  of  the  domain.  The  nodes  are  free  to  move  vertically  by 
symmetric  boundary  conditions  along  the  left  and  right  edges  of  the  domain. 

This  problem  can  be  solved  readily  by  perhaps  many  PDE  solution  methods 
whenever  v  assumes  sufficiently  large  values.  For  example,  a  value  of  v  » 
0.02  produces  shock  gradients  on  the  order  of  10*. 

The  MFE  method  requires  only  an  8  x  8  grid  to  give  reasonably  accurate 
solutions  to  this  problem,  and  Figures  3  and  4  show  very  accurate  MFE  solu¬ 
tions  on  a  12  x  12  grid.  Here  Figure  3  presents  an  Isometric  view  of  the 


I  i  i  i  i  i  i  i  i  s  i 

X 

Figure  5.  Contour  plots  of  selected  values  of  u  at  t  =  3.0  in  the  2-D 
Burger-like  example  on  a  12  x  12  grid. 


evolving  profile  of  the  y  component  of  velocity  at  t  *  3.0,  well  after  the 
shock  has  formed  and  after  the  wavefront  has  undergone  significant  shearing. 
The  x-component  of  velocity  Is  sufficiently  sheared  that  a  hidden  line  plot, 
which  Is  not  yet  available.  Is  required  for  easy  Interpretation  by  the  naked 
eye.  The  MFE  grid  nodes  have  migrated  extensively  from  their  Initial  posi¬ 
tions  as  can  be  seen  In  Figure  4  which  represents  the  grid  mesh  projected 
onto  the  x-y  plane  at  t  *  3.0.  Figures  5  and  6  present  contour  plots  for 
selected  constant  values  of  u  and  v,  respectively,  at  t  *  3.0.  It  Is  evident 
from  the  magnitudes  of  shock  gradients  and  from  the  regions  of  significant 
curvature  which  span  nearly  the  entire  domain  that  an  alternative  POE  method 
with  a  fixed  grid  may  require  on  the  order  of  10*,  or  more,  grid  nodes  In 
order  to  achieve  comparable  degrees  of  accuracy  In  this  problem. 

This  same  basic  problem  can  now  be  made  to  correspond  to  a  much  more 
demanding  physical  problem  by  letting  v  *  0.002.  Figure  7  shows  an  Isometric 
view  of  the  MFE  solution  on  a  16  x  16  grid  for  this  case.  Shock  gradients 
are  now  generated  with  magnitudes  of  several  times  103.  Before  discussing 
these  MFE  results  In  detail,  some  general  observations  should  be  discussed: 
It  Is  extremely  unlikely  that  any  other  existing  POE  method  using  either  a 
fixed  grid  or  a  less  than  optimal  adaptive  grid  can  accurately  solve  this 
test  problem  with  fewer  than  10®-10®  grid  nodes.  It  should  also  be  noted  here 
that  nunerous  Invlscld  solvers  which  are  under  development  do  not  apply  at 
all  to  this  type  of  advect Ion-diffusion  problem  because  the  Laplacian  Is  an 
essential  mathematical  operator  whose  effects  must  be  rigorously  resolved  In 
advectlon-dlffuslon  PDE's.  Because  Invlscld  solvers  do  not  generally  solve 
PDE's  which  contain  Laplaclans,  they  generate  shocks  with  gradient  shapes  and 
magnitudes  that  are  governed  exclusively  by  the  selected  grlddlng  and/or  by 
the  purely  numerical  dissipative  processes  In  the  Invlscld  method,  per  se. 
Consequently,  Invlscld  solvers  have  no  chance  of  resolving  correctly  any  of 
those  Physical  dissipation  effects  which  are  usually  expressed  by  Laplacian 
operators  and  are  present  with  fundamental  physical  significance  In  transport 
theory,  hydrodynamics,  plasma  physics,  continuum  mechanics,  and  many  other 
disciplines  In  the  physical  scler^es.  Thl  critical  discussion  Is  not  In¬ 
tended  to  denigrate  the  extensive  <?*ea  n  efforts  on  Invlscld  POE  solvers 
and/or  fixed  node  POE  methods;  but  it  does  suggest  that  efforts  to  accommo¬ 
date  Laplacian  operators  In  otherwise  Invlscld  solution  methods  and  efforts 


Figure  7.  Isometric  view  of  Burger's  test  example  at  t  «  1.8  with 

v  ■  0.002  on  a  16  x  16  MFE  mesh.  Shock  gradients  have 

magnitudes  of  approximately  103  In  this  solution  for  u. 

to  Investigate  more  optimal  adaptive  grid  methods  for  use  In  many  existing 
PDE  methods  which  are  applied  to  advectlon-dlffuslon  problems  should  now 
assune  greatly  Increased  significance.  In  the  meantime,  the  WE  method  Is 

proving  to  be  a  certain  kind  of  research  pacesetter,  and  It  Is  providing 

various  clues  to  some  of  the  significant  new  areas  ihere  mathematics  research 
can  profitably  be  Intensified,  as  will  be  discussed  further  In  other  task 
areas  below. 


Early  MFE  results  In  2-D  are  apparently  continuing  the  trend  which 
appeared  In  previous  1-0  results.  There,  MFE  solutions  of  both  the  Navler- 
Stokes  and  physically  dissipative  continuum  mechanics  equations  In  1-D  ex¬ 
hibited  perhaps  unprecedented  simultaneous  resolution  of  extremely  disparate 
microscale  and^  macroscale  physical  processes. While  2-0  MFE  results 
which  emerged  during  this  period  exhibit  similar  promising  features,  numerous 
mathematical  problems  still  require  resolution  In  order  to  attain  fully  the 
desired  levels  of  success  in  truly  large-scale  problems  in  2-D.  Clues  to 
these  problem  areas  can  be  seen  In  Figure  7.*  For  example,  the  Irregularity 
of  the  grid  triangles  In  the  face  of  the  shock  could  eventually  prove  to  be 
troublesome.  Similarly,  the  small  oscillation  at  the  base  of  the  shock  In 
this  run  Is  unsatisfactory,  even  though  it  can  be  eliminated  in  any  number  of 
ways.  Extensive  testing  and  analysis  has  indicated  that  the  causal  mechanisms 
underlying  such  mesh  Irregularities  and  oscillations  In  2-D  can  be  associated 
with:  (1)  time  step  and  error  control  policies  in  the  basic  ODE  integrator 
of  Gear  which  Is  presently  used,  (11)  convergence  properties  of  the  linear 
solver,  and  (111)  limitations  In  the  first-generation  regularization  func¬ 
tions.  Each  of  these  areas  has  been  under  Intensive  investigation  during  this 
period  and  some  early  results  and  their  Implications  are  discussed  below. 

2.  ODE  SOLVERS  FOR  POE  METHODS 

The  current  status  In  this  task  area  Is  that  most  existing  ODE  solvers 
are  not  well-suited  for  ready  Implementation  In  either  the  MFE  method  or 
nunerous  other  advanced  PDE  methods.  This  critical  comment  Is,  again,  not 
Intended  to  denigrate  the  Impressive  advances  In  ODE  research  and  development 
during  the  past  decade;  Instead,  It  Is  Intended  to  bring  a  strong  new  focus 
upon  the  needs  of  PDE  solution  methods.  In  general,  and  more  specifically 
upon  the  pressing  needs  of  adaptive  grid  PDE  methods  which  may  Involve  large 
numbers  of  discretized  equations  with  highly  distorted  grids.  Large  distorted 
grid  meshes  may.  In  turn,  augur  for  Iterative  linear  solvers  which  can  solve 
poorly  conditioned  matrix  equations,  as  will  be  discussed  further  below. 


*It  Is  apparent  that  these  suggested  mathematical  problems  will  have  to  be 
resolved  not  only  for  the  MFE  method  but  also  for  most  other  advanced  PDE 
methods  which  may  seek  to  solve  the  difficult  advect Ion-diffusion  equations 
which  frequently  arise  In  physical  problems. 
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A  basic  difficulty  with  most  stiff  ODE  software  has  come  to  the  fore  In 
the  present  MFE  research;  l.e.,  most  ODE  solver  packages  have  been  designed 
to  accommodate  many  different  types  of  classic  ODE  problems.  By  classic,  one 
refers  here  to  such  problems  as  chemical  kinetics  systems  In  which  the  depen¬ 
dent  variables  (e.g.  species  concentrations)  are  all  generlcally  the  same. 
The  error  and  time  step  controlling  policies  In  solvers  for  classic  ODE  pro¬ 
blems  are  usually  less  than  satisfactory  for  applications  of  the  ODE  package 
to  PDE  solution  methods.  In  PDE  systems,  the  spatial  dependence  of  generl¬ 
cally  dissimilar  variables  comes  Into  play.  In  fluid  dynamics  problems,  for 
example,  the  overall  array  of  PDE  variables  which  have  been  discretized  on  N 
grid  nodes  (x^,  x2,  ...  xN)  can  be  represented  as  {pj,  m^,  Ej,  p2,  m2,  E2,... 

PlfN*  hh  where  pj  =  pfxj),  p2  =  p(x2),  etc.  An  ODE  solver  then  operates  on 
this  array  of  discretized  PDE  variables  as  a  single  large  vector  (y^  y2,  y3, 

y4,  ...  y3N_2-  *3IM'  V,"l'ere  yi  ■  Pi.  *  «i.  y3  -  Ei.  y4  •  p2- etc- 

cause  the  error  control  policies  In  the  Gear  ODE  package,  for  example,  are 
based  upon  L2  norm  of  all  normalized  quantities  y^/(y^ )max  unacceptably  large 
errors  can  be  admitted  In  some  Individual  components  of  p,  m,  or  E  at  arbi¬ 
trary  spatial  locations.  A  much  better  measure  for  error  control  policies  in 
PDE  applications  are  maximum  norms  applied  to  each  discretized  PDE  variable. 
The  Implementation  of  alternative  norms  Is  found  to  extend  deeply  into  the 
logical  structure  of  most  ODE  software  packages,  and  alterations  must  usually 
be  performed  by  someone  who  Is  Intimately  familiar  with  the  ODE  package. 
(Dr.  Said  Doss  Is  the  local  expert  on  these  matters  in  our  MFE  research.) 

In  view  of  such  considerations,  significant  levels  of  effort  have  been 
devoted  during  this  work  period  to:  (1)  modifications  of  Gear's  basic  ODE 
method  for  MFE  computations;  this  Involved  wholesale  alterations  of  the 
Internal  Gear  code  structure  and  also  extensive  considerations  of  scaling  of 
MFE  problem  variables;  and  (11)  Inevitably,  the  development  of  entirely  new 
ODE  Integraton  procedures  which  better  serve  PDE  solution  needs. 

The  extensive  modifications  to  Gear's  method  have  sufficed  to  solve 
moderately  challenging  POE's  with  modest  nunbers  of  MFE  grid  nodes  as  was 
seen  above.  But  It  Is  now  clear,  also,  that  completely  new  ODE  code  struc¬ 
tures  will  be  needed  In  pending  large-scale  MFE  computations.  We  have, 
therefore,  undertaken  the  development  of  a  low-order  Runge-Kutta  Integration 
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package  for  MFE  computations.  This  solver  addresses  several  PDE  needs:  First, 
error  control  measures  operate  on  flexibly  ordered  variable  arrays  using 
maximum  norms  on  a  (POE)  component-by-component  basis.  Second,  PDE  solutions 
have  been  found  to  require  much  more  gradual  time  step  advancement  policies 
than  have  been  built  Into  most  classic  ODE  solvers.*  This  distinction  between 
classic  ODE  and  PDE  time  step  properties  apparently  stems  from  the  fundamen¬ 
tally  coupled  space-time  dependences  In  PDE  systems,  vis  a  vis  classic  ODE 
problems  which  have  no  direct  or  Implied  spatial  dependences.  Whereas  It  Is 
computationally  worth  the  effort  to  attempt  very  large  Incremental  Increases 
(sometimes  by  several  orders  of  magnitude)  In  At  In  classic  ODE  applications 
—even  If  such  attempts  may  sometimes  fall --one  finds  that  the  computational 
penalties  for  unsuccessful  large  At  Increases  are  much  more  severe  In  those 
PDE  applications  where  space-time  couplings  augur  Intrinsically  for  more 
gradual  At  advancement  policies.  Third,  time  step  control  policies  In  the 
new  ODE  solver  also  Incorporate  convergence  criteria  from  Iterative  linear 
systems  solvers.  Such  Iterative  solvers  should  henceforth  be  used  In  large- 
scale  MFE  computations  In  the  Interest  of  minimizing  computer  memory  require¬ 
ments.  Fourth,  low-order  ODE  methods  are  now  used  because  high-order  solvers 
provide  no  apparent  advantages  In  MFE  applications  and  because  low-order 
methods  simplify  the  numerical  logic.  Improve  the  code  reliability,  and  avert 
possible  errors  associated  with  changes  of  order  which  are  sometimes  present 
In  classic  ODE  system  solvers.  Finally,  constraints  on  allowable  fractional 
changes  In  PDE  dependent  variables  are  Incorporated  In  the  overall  time  step 
control  policy  In  the  new  ODE  solver.  This  new  ODE  solver  Is  presently  being 
Implemented  for  use  In  large-scale  MFE  computations;  and  detailed  descriptions 
of  this  solver.  In  conjunction  with  MFE  test  applications,  will  appear  In 
forthcoming  journal  submissions. 

From  these  Initial  results  It  Is  clear  that  renewed  ODE  research  efforts 
on  PDE-related  problems,  from  several  possible  conceptual  bases.  Is  now  time¬ 
ly  If  not  long  overdue.  Our  own  efforts  have  barely  begun  to  uncover  many  of 
the  most  pressing  needs,  much  less  to  perform  the  extensive  detailed  tasks  of 
nunerlcal  analysis  which  should  now  be  pursued.  Certainly,  the  present  MFE 
project  would  benefit  from  expanded  ODE  efforts  which  deal  with:  (1)  scaled 


*These  policies  also  extend  deeply  Into  the  ODE  code  structure. 
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and  unsealed  systems  of  Incommensurate  variables  on  arbitrarily  connected 
grid  meshes*  (11)  linkage  of  certain  At-sensltlve  convergence  criteria  Into 
the  general  ODE  Integration  control  policy*  and  (111)  splitting  of  the  solu¬ 
tion  of  the  grid  node  equations  from  the  solution  of  the  discretized  equa¬ 
tions  for  the  physical  variables. 

3.  LINEAR  SOLVERS  FOR  THE  MFE  METHOD 

Advect Ion-diffusion  equations  have  steadfastly  resisted  (If  not  defied) 
satisfactory  numerical  solution  whenever  they  have  been  used  to  describe 
physical  processes  over  highly  disparate  scales.  Such  problems  occur,  for 
example.  In  numerous  applications  of  Navler-Stokes  equations  to  viscous  com¬ 
pressible  fluids  which  may  contain  shear  layers,  shocks,  and  separated  flows. 
The  basic  difficulty  derives  from  the  nature  of  the  matrix  equations  which 
must  be  solved  In  numerical  PDE  methods  that  are  applied  to  these  problems. 
The  matrix  equations  for  discretized  advectlon-dlffuslon  PDE's  are  large, 
sparse  linear  systems  In  which  the  matrices  are  non-symmetrlc  and  are  not 
dominated  by  terms  on  the  diagonal.  The  skewness  of  these  PDE  matrices  can 
become  quite  large  for  large  At's  and  for  highly  distorted  grid  meshes,  both 
of  which  are  key  factors  In  efficient  solutions  of  these  types  of  advectlon- 
dlffuslon  equations. 

The  simple  advectlon-dlffuslon  equation,  yt  +  £  •  Vy  *  W^y,  can  be  used 
to  Illustrate  some  of  these  features.  Upon  discretization  of  the  advectlon- 
dlffuslon  equation,  a  linear  system  of  the  form  (A+B)X  *  C  must  be  solved. 
The  matrix  A  represents  the  stencil  associated  with  the  advectlon  operator, 
and  the  matrix  B  represents  the  stencil  associated  with  a  nine  point  differ¬ 
ence  scheme  for  the  Laplaclan  operator  In  2-D.  For  representative  values  of 
At,  these  matrices  may  contain  elements  with  the  scaled  magnitudes  shown 
below: 
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Testing  and  analysis  has  revealed  that  most  available  linear  solvers  have 
^  relatively  poor  rates  of  convergence  when  such  significant  large  elements  can 

P  occur  away  from  the  diagonal  In  non-symmetrlc  matrices.  For  example,  such 

Iterative  matrix  solution  methods  as  conjugate  gradient,  multi-grid  and 
numerous  other  modern  linear  systems  solvers  which  work  well  for  symmetric 
matrices  In  discretized  elliptic  equations  and/or  for  uniform  grid  meshes  do 
P  not  converge  satisfactorily  In  presently  considered  advect Ion-diffusion 

problems.  Hie  source  of  difficulty  for  the  existing  linear  solvers  clearly 
derives  from  the  highly  skewed  matrices  and  their  off-diagonal  dominance.  We 
have  also  shown  during  this  work  period  that  the  direct  L-U  decomposition 
7?  method  which  has  been  used  In  the  Gear  method  until  recently  becomes  both 

noisy  and  computer  storage  limited  when  large  bandwldths  arise  In  problems 
with  more  than  a  moderate  number  of  MFE  grid  nodes.  Again,  we  do  not  wish  to 
m  denigrate  the  extensive  ongoing  work  on  linear  systems  solvers— but,  rather, 

to  call  attention  to  these  essential  keys  to  further  progress  in  certain  PDE 
research  areas. 
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Having  Identified  more  clearly  the  significance  of  these  issues,  we  have 
developed  In  collaboration  with  Professor  Keith  Miller  one  promising  new 
approach  to  handling  more  effectively  these  Imposing  PDE  requirements  on 
linear  systems  solvers.  This  new  matrix  solution  scheme  has,  so  far,  achieved 
good  convergence  rates  for  At's  which  may  be  10  to  20  times  greater  than  the 
large  values  of  At  called  for  by  the  ODE  Integrator.  (It  Is  generally  hoped 
In  PDE  solutions  that  the  time  step  size  Is  determined  by  the  truncation 
error  of  the  ODE  Integrator  and  not  by  severely  limited  convergence  proper¬ 
ties  of  the  linear  solver.)  This  advanced  linear  solver  has  solved  the 
Burger's  equations  discussed  above  with  the  same  CPU  cost  as  the  direct  L-U 
decomposition  method  In  the  Gear  solver,  but  with  greatly  reduced  storage 
requirements.  Implementation  of  this  new  linear  solver  for  large-scale  MFE 
computations  Is  progressing  well,  and  details  of  this  new  approach  will 
appear  In  journal  submissions  during  the  next  work  period.  Clearly,  our 
Initial  efforts  on  a  more  adequate  linear  solver  for  advectlon-dlffuslon 
problems  have  only  barely  opened  a  new  facet  of  research  which  now  warrants 
much  more  Intensive  theoretical  and  practical  analysis,  both  by  ourselves  and 
by  many  others. 
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4.  REGULARIZATION 


Regularization  techniques  have  rarely  been  used  systematically.  If  at 
all.  In  PM  research  In  the  past.  There  Is  thus  presently  great  confusion 
and  misunderstanding  of  the  role  of  regularization  techniques  In  POE  solution 
methods.  On  the  one  hand,  many  practitioners  of  conventional  PDE  methods 
suspect  that  regularization  Is  an  unfair  trick  by  which  one  can  force  PDE 
solutions  to  come  out  In  any  desired  manner  and  that  such  methods  can  there¬ 
fore  not  be  trusted.  On  the  other  hand,  regularization  methods  are  proving 
to  be  valid  and  powerful  mathematical  tools  which  can  now  be  applied  to 
achieve  effective  grid  movement  criteria  systematically  and  to  ensure  that 
high  PDE  solution  accuracy  Is  also  achieved  In  the  process. 

Only  the  simplest,  first-generation  regularization  functions  have  been 
used  In  2-D  MFE  problems  to  date.  These  penalty  functions  act  like  springs 
and/or  dashpots  In  their  action  on  mesh  triangle  altitudes  (see  Figure  8). 
The  current  penalty  functions  allow  mesh  triangles  to  distort  more  or  less 
arbitrarily,  so  long  as  altitude  magnitudes  remain  positive  and  maintain  some 
designated  minimum  separation.  This  simple  strategy  has  worked  remarkably 
well  In  a  surprisingly  broad  range  of  2-D  problems.  Including  all  of  the 
results  shown  In  this  report.  We  are  presently  probing  the  limits  of  ade¬ 
quacy  of  this  simple  first-generation  regularization  method  In  such  extended 
applications  as  Mach  stems  In  gas  dynamics  and  dynamic  impact  problems  In 
continuum  mechanics.  In  the  Burger's  test  example  discussed  above  with  v  * 
0.002  and  a  16  x  16  MFE  mesh,  minimum  nodal  separations  were  chosen  to  be 
much  smaller  thanv.  As  seen  In  earlier  figures,  the  grid  triangles  can 
become  very  Irregular  and  assume  configurations  with  extremely  large  aspect 
ratios  (O(102)  to  O(103)).  Figures  9  and  10  below  show  this  same  Burger's 
test  problem  run  with  "softer"  dashpot  forces  acting  on  triangle  altitudes 
than  In  the  run  shown  previously  In  Figure  7.  The  triangles  In  the  face  of 
the  shock  are  more  regular  In  this  latter  run,  and  their  compaction  In  the 
center  of  the  wavefront  resolves  the  region  where  Intense  shearing  occurs. 
Figure  11  shows  the  projection  of  the  MFE  grid  mesh  on  the  x-y  plane  In  this 
example.  Selected  mesh  connections  have  been  traced  In  heavy  Ink  In  these 
latter  figures  In  order  to  Indicate  the  general  migration  pattern  of  the  MFE 
nodes.  It  Is  evident  that  extremely  large  mesh  triangle  aspect  ratios  have 
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Figure  8.  Schematic  Representation  of  First -Generation  MFE 
Regularization  Functions  In  2-D. 

In  the  long  run,  there  Is  no  fundamental  reason,  or  desire,  for  the  MFE 
mesh  to  always  be  as  highly  skewed  as  the  physical  flow  lines  In  order  to 
accurately  resolve  such  shear  layers.  (It  Is  nevertheless  encouraging  at 
this  stage  of  development  that  the  mathematical  potency  of  the  MFE  method  Is 
sufficiently  great  to  handle  such  large  grid  aspect  ratios  effectively.)  Al¬ 
ternative  regularization  criteria  which  would  promote  mesh  homogeneity  (e.g., 
by  minimization  of  grid  triangle  aspect  ratios)  are  presently  under  consider¬ 
ation  for  use  In  conjunction  with  the  first-generation  regularization  func¬ 
tions.  Several  possible  formulations  of  grid  homogenizing  regularization 
schemes  will  be  tested  In  the  next  work  period.  These  new  criteria  are  ex¬ 
pected  to  improve  nunerlcal  conditioning  properties  and  thus  lead  to  greater 
computational  economy  In  MFE  codes.  Mesh  homogenlety  should  also  be  a  signi¬ 
ficant  factor  In  applications  which  must  resolve  turbulent  eddies  and/or 
other  rotating  flows. 
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Figure  9.  Isometric  view  of  the  2-0  solution  of  the  velocity 

component,  v.  In  the  Burger-like  example  with  v  ■  0.002 
on  a  16  x  16  NFE  grid. 
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Figure  10.  Isometric  view  of  the  2-D  solution  of  the  velocity 

component,  u,  In  the  Burger-like  example  with  vB  0.00 
on  a  16  x  16  MFE  grid.  Note  that  the  viewing  angle  Is 
rotated  by  90*  for  a  clearer  view  of  the  doubly  skewed 
wavefront . 
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We  have  now  barely  opened  a  potentially  vast  area  of  Investigation  of  re¬ 
gularization  techniques  for  adaptive  POE  meshes.  This  is  certainly  one  of  the 
key  areas  where  we,  and  perhaps  many  others,  can  profitably  expend  greater 
efforts,  particularly  In  view  of  the  very  attractive  MFE  properties  which  are 
emerging  In:  (1)  alternative  co-ordinate  systems  and  (11)  the  treatment  of 
Interface  phenomena,  both  of  which  have  a  direct  bearing  on,  and  relationship 
to  alternative  regularization  methods.  Here  also  the  MFE  method  Is  a  certain 
kind  of  research  pacesetter  which  Is  suggesting  means  of  applying  these  new 
mathematical  methods  of  regularizing  POE  node  motions  In  both  new  and  conven¬ 
tional  PDE  methods  so  as  to  resolve  such  historically  persistent  dilemmas  as 
singularities  at  origins  of  spherical  and  cylindrical  coordinate  systems, 
artificial  smearing  of  Interfaces,  unduly  constrained  grid  aspect  ratios, 
numerical  diffusion,  severe  time  step  constraints,  and  grid  node  utilization. 

5.  ALTERNATIVE  CO-ORDINATE  SYSTEMS 

Initial  work  Is  In  progress  on  2-D  MFE  calculations  In  cylindrical  co¬ 
ordinates.  The  results  In  this  task  also  provide  guidance  for  later  develop¬ 
ments  In  spherical  co-ordinates.  The  major  issue  of  present  interest  Is  the 
apparently  natural  elimination  of  singularities  at  the  origin.  Such  singu¬ 
larities  have  historically  plagued  many  conventional  PDE  methods. 

Transport  equations  contain.  In  cylindrical  co-ordinates,  for  example, 

1  3 

advectlon  operators  of  the  form  (ry),  where  r  Is  the  radial  co-ordinate; 
and  y  Is  a  dependent  variable  of  the  PDE  system.  Singularities  or  other 
anomalous  features  frequently  occur  In  various  discretized  representations  of 
the  term  (y/r)  as  r— 0.  In  contrast,  the  MFE  discretization  Is  formulated 
In  terms  of  well-defined  Inner  products  which  eliminate  such  possible  singu¬ 
larities.  For  example,  the  Inner  products  of  the  term  (y/r)  with  the  basis 
function  a,  taken  over  the  Interval  Ar,  Is  given  by 

J  (y/r)  «  a  •  rdr  ■  I  y  •  otdr  . 

A  r  A  r 

The  integral  of  a  •  y  dr  Is  essentially  analytic  and  Is  readily  evaluated 
everywhere  on  the  problem  domain.  This  attractive  MFE  property  In  cylindrical 
co-ordinates  obviously  holds  In  a  similar  manner  In  spherical  co-ordinates. 
The  properties  of  these  r-welghted  norms  are  naturally  different  than  the  MFE 


Inner  products  which  Mere  used  In  the  Cartesian  co-ordinate  systems  considered 
in  earlier  MFE  work.  Analysis  and  testing  of  these  properties  associated  with 
r-welghted  norms  and  of  node  controlling  penalty  functions  In  cylindrical  co¬ 
ordinates  have  been  undertaken  In  this  period,  and  extensive  work  Is  expected 
to  continue  In  this  area  In  the  future  In  order  to  understand  and  exploit 
fully  the  benefits  of  this  analytic  MFE  formulation  of  otherwise  troublesome 
PDE  operators  In  cylindrical  and  spherical  co-ordinates. 
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SEMINARS  GIVEN  ON  OUR  MFE  RESEARCH 

1.  Air  Force  Weapons  Laboratory  (Major  Rajmond  Bell) 

2.  Sandla  Livermore  Laboratory  (Dr.  R.  Kee) 

3.  University  of  Toronto  (Prof.  I.  Glass) 

4.  Lawrence  Livermore  Laboratory  (H.  Division  Seminar) 

5.  Exxon  Research  (Dr.  G.  Byrne) 


CONSULTIVE  FUNCTIONS  TO  OTHER  AGENCIES 

1.  Eglln  Air  Force  Base;  we  are  currently  under  contract  to  Investigate 
the  application  of  the  WE  method  In  2-D  to  armor  penetration  effects;  con¬ 
tact:  Major  Guy  Spltale. 

2.  Defense  Nuclear  Agency;  we  are  currently  under  contract  to  Investi¬ 
gate  the  application  of  the  MFE  method  In  2-D  to  alrblast  effects;  contact: 
Dr.  George  Ullrich. 


